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. Abstract. The impurity dynamics in stellarators has become an issue of moderate 

' concern due to the, a priori, inherent tendency of the impurities to accumulate in the 

core when the neoclassical ambipolar radial electric field points radially inwards (ion 
root regime). This accumulation can lead to collapse of the plasma due to radiative 
</} 1 losses, and thus limit high performance plasma discharges in non-axisymmetric devices. 

Theoretically, a quantitative description of the neoclassical impurity transport is 
&\ complicated by the breakdown of the assumption of small q&/T for impurities, where 

q is the electric charge, T the temperature in energy units, and 4> the electrostatic 
O ( potential variation within the flux surface. 

The present work describes quantitatively the particle transport of impurities in the 
■ frame of local neoclassical theory when q<t/T = 0(1) in the Large Helical Device (LHD) 

^ \ stellarator. The central numerical tool used is the Sf particle in cell (PIC) Monte Carlo 

£^ ■ code EUTERPE. The 4> used in the calculations is provided by the neoclassical code 

. GSRAKE. The possibility of obtaining a more general self-consistently with EUTERPE 

t-H | is also addressed and a preliminary calculation is presented. 

OV 

o: 

T"! , 1. Introduction 

>■ 

• i-H ■ 

^ ; Thermonuclear fusion would benefit from the achievement of quasi steady-state magnetic 

e$ ; plasmas confinement with similar characteristics to those expected in a future reactor. 

In this respect the stellarator concept has an advantage over the pulsed tokamak. On the 
other hand the neoclassical transport exhibited by the former at low collisionality in the 
absence of electric fields, in the so-called \ jv regime, is considerably larger. The reason 
for this unfavourable behaviour is, that in contrast to the axisymmetric case, where 
the collisionless trajectories of the trapped particle are confined, the orbits of particles 
trapped in the helical magnetic wells are generally not confined. This situation leads to 
a different radial transport rate, for each species in the plasma, and necessitates a radial 
electric field E r = — V$o that restores ambipolarity. $ is the lowest order electrostatic 
potential and only depends on the flux surface label. In the present work we use the 
normalised toroidal flux s = r^ ff , r e ff = r/a the effective radius and a the minor radius of 
the plasma. In the present context the basic ordering parameter, S = p/L, is the Larmor 
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radius p normalised to a typical macroscopic variation length scale L. Consequently the 
distribution function is expressible up to first order as / « f + 5/, with Sf/fo ~ 0(5). 
The sign of E r = |E r | in standard conditions is negative (ion root regime) and predicts 
accumulation of impurities pQ. The subsequent increase of radiative losses from the 
core can cause the collapse of the plasma (2j [3], and, in the worst case, endanger the 
capability of the stellarator to confine it in steady-state. 

A quantitative and comprehensive description of the impurity dynamics requires not 
only the ambipolar radial electric field but also the electrostatic potential $ determined 
by 5f. The usual neglect of this portion of electrostatic potential for bulk ions and 
electrons rests on two assumptions. First, the radial E x B drift arising from $, 

V$xb 

v 4> = J} , (!) 

is assumed to be smaller than the curvature and grad-B drifts. At low magnetic pressure 
(/?—»► 0) these are given by, 

v d = ~—B2 Ah x V5 > ( 2 ) 
q B z 

with m the mass of the particle, q = Ze its electric charge, Z the charge number, e the 
unit charge absolute value, B the magnetic field strength, b a unitary vector pointing 
in the direction of the magnetic field line, /i = v\/2B the normalised magnetic moment, 
v\\ the parallel velocity and v± the perpendicular one. It can be shown that the ratio 
between the absolute values of v$ and is thus of order 

v# _ Ze$ R ^ 



v& T a 

where R is the major radius of the device, and it is assumed that the typical variation 
length scales for B and $ are similar to R and a, respectively. This ratio can be 
acceptably small at low values of Z but is usually considerable for heavy impurities. 
Second, the parallel acceleration a\\ = —(q/m)b • V$ is assumed to be negligible 
compared to the mirror force a m = — fib • VB. The ratio between these two is of 
order, 

a\\ Ze§ B 

^ ~ ~T^AB' ^ ' 

with AB the typical amplitude of the helical magnetic wells. Again, the proportionality 
to Z makes necessary to account for a\\ if the impurity transport is to be quantitatively 
computed, although $ can be sufficiently small for the bulk species. 
On the other hand, the inclusion of $ into the problem makes the kinetic energy of 
the impurities vary enough to violate the traditional neoclassical theory mono-energetic 
modus operandi. 

In the present work we put the focus on the computation of the neoclassical particle flux 
of impurities including $, thus abandoning the mono-energetic assumption, but keeping 
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the radially local one. The calculations are performed with the Monte Carlo 5f PIC code 
EUTERPE. A concise description of it is given in section [2J highlighting the truncation 
of the global characteristics to perform local neoclassical runs. Section [3] shows particle 
flux calculations for C 6+ , Ne 8+ and Fe 20+ in the LHD standard configuration including 
the poloidal variation of with 6 the poloidal coordinate. This is obtained from 

the solution of the steady-state ripple averaged drift kinetic equation obtained by the 
GSRAKE code [3j. Section [4] addresses the calculation of $ by EUTERPE, discussing a 
preliminary result that includes the dependence of $ on the toroidal coordinate (f) and 
time t. Finally, in section [5] a summary of the results and a discussion are presented. 



2. The EUTERPE code in the local neoclassical limit 

EUTERPE is a global 5f PIC Monte Carlo code, full-radius and full-flux surface, 
initially conceived for numerical simulation of linear gyro-kinetic micro-turbulence in 
3D equilibria (5j El [7] . After undergoing successive updates the current version is non- 
linear, can treat multiple kinetic species simultaneously, and perform electro-magnetic 
simulations. Recently it has been extended to include pitch angle scattering collisions 

® 

The set of phase space coordinates that EUTERPE uses is z = {R, v\\ , /i}. R is the position 
of the guiding center of the particle in neoclassical runs, or its gyro-center in gyro- kinetic 
ones. The collisionless trajectory of a particle in phase space in the electrostatic limit, 
for simplicity written at /3 — 0, is determined by the following set of equations: 



R = v,| +V£ + v d , (5) 

a v\\ 
— b • E — /ib • VB + — | 

m B z 

ft = 0, (7) 



— b • E — /xb • VB + — !j- (b x VB) ■ E, (6) 



where an overhead dot denotes a time derivative, E = — V$ is the electric field, 
v^ = Exb/B and $ the electrostatic potential obtained from the gyro- kinetic quasi- 
neutrality equation with arbitrary spatial dependence. The numerical splitting of 
the distribution function assumes that the lowest order part is a local Maxwellian, 

fo = f M (s,v\\,v±) = {2n ) n sil\ h ( s ) ex P > with v th = V T / m and T = T ( s ) 

the temperature. The departure from it, 5/(z,t), is simulated using markers, whose 
trajectories are pushed according to eqs. ©-(CD), periodically interrupted by the 
application of a random change of their pitch angle to account for collisions [9]. 
It is straightforward to prove that eqs. ©-(CD) conserve the total energy of the 
particle and preserve the incompressibility of the phase space flow, i.e. 8 = with 
8 = v\j2 + jiB + (q/m)$ and V z • z = with V z = (V, d V{{ , . 

The trajectory in real space followed by a particle given by eq. ([5]) is referred to as 
global since the radial magnetic and E x B drifts across the flux surface, • Vs and 
v E - Vs respectively, are accounted for. In contrast, the guiding center trajectory without 
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Figure 1: Particle flux of hydrogen in an axisymmetric circular cross section tokamak for 
the radial density flat profile of n = 10 19 m -3 and a temperature profile T — 2(1 — r| ff ) 
keV (left), and for n = (1 — rf ff ) x 10 19 m -3 at a fixed temperature of T = 2 keV (center). 
In the right, particle flux density for C 6+ in LHD for the profiles n = (5 — 0.2r^) x 10 19 
T = (3 — 0.93rg ff ) keV and the E r profile embedded in the plot. 



these drifts in lowest order is called local. Local trajectories find routine application in 
neoclassical simulations and follows from assuming E ~ E r , resulting in laying on the 
flux surface, and v&/v\\ ~ S. For impurities we need to retain $ though, i.e E = E r — V$, 
and assume that v$/v\\ ~ S too. In this approximation energy is not conserved since 
the neglect of unbalances a cancellation in £ resulting in £ — qv^ • E r . The variation 
of the energy during one collisional time r can be estimated to be lower than a certain 
value, (A£) T /T < (g$ /T)/(Ar/a), with Ar < v&t the scale of the total radial drift in 
a typical collisional time, and V ^ a -1 . 

Our equations for the marker trajectory in phase space for a local neoclassical run are 
thus, 



R = V||+V£, (8) 

m = -^b • V$ - fib • VB + (b x VB) • E r , (9) 

rn B z 

A = o, (10) 

with = E r x h/B. For the tokamak, where the small- E r limit ve/vu ~ S applies pU] 
(unless the plasma rotates rapidly) the second term in eq. (JSj) and the third in eq. (|9j) 
are neglected. The resulting kinetic equation finally reads as follows: 



w + ^ . w/ + .m = _ (va + } . v/m+ 

at ov\\ , . 

, (11) 

+ -#(v||+v d )-(E r -V$). 

mv z h v 11 7 V / 

EUTERPE, in the present neoclassical modality, was satisfactorily benchmarked against 
theory and the codes DKES and GSRAKE. Figured] shows a comparison between the particle 
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Figure 2: (a) Ambipolar radial electric field obtained with GSRAKE in the two cases 
considered for the impurity runs launched by EUTERPE. The corresponding $ are 
represented in plot (b) for case I, and (c) for case II. 



fluxes obtained with these codes in different cases. It is important to remark that in this 
comparison $ is missing since DKES and GSRAKE assume mono-energetic trajectories. 

3. Impurity particle transport in the presence of $ (9) 

The code GSRAKE solves the ripple- averaged drift kinetic equation providing, apart from 
the neoclassical fluxes, the first order corrections to the equilibrium density and elec- 
trostatic potential, h (r, 9) and $ (r, 9). Since the code assumes mono-energeticity the 
impurity fluxes are not obtained. The magnetic configuration is accounted for by means 
of a multiple-helicity model, and the bounce average is performed along the toroidal co- 
ordinate (j) making the dependence on it vanish. The resulting ^-dependent n and $ are 
written in Fourier series, and their coefficients are iteratively adjusted to simultaneously 
fulfil the quasi-neutrality and ambipolarity constraints among bulk ions and electrons 
(for the details of the code see refs. [U [TT]). 

In the present section the E r and <E> obtained with GSRAKE for two different sets of density 
and temperature profiles have been used as input for EUTERPE, which in turn obtains 
the particle flux F z . The magnetic equilibrium considered corresponds to the standard 
LHD configuration, R « 3.73 m, a « 0.6 m and £>oo(?~eff = 0.5) = 2.53 T. Two cases 
have been considered and are labelled as case I and case II The equilibrium tempera- 
ture and density profiles considered in the former case for the bulk ions (hydrogen) are 
n = (6 — 5r^ ff ) x 10 19 m -3 and T = 2(1 — 2.8r^) keV. The rather lower collisional case 
II considers instead n = (5 — 0.2r|! ff ) x 10 19 and T = (3 — 0.93rg ff ). The temperature 
profiles are taken to be the same for electrons and impurities, and their central density 
is adjusted to fulfil quasi-neutrality. The effective charge is set to — 1.05 at all radii. 
In fig. [2](a) the ambipolar radial electric field is displayed for both pairs of profiles, in 



Neoclassical impurity transport in stellarator geometry 



0.0 
-1.0 
^ -2.0 
1 -3.0 

CL 

\ -4.0 
-5.0 
-6.0 
-7.0 







W/O <I> 

with $ 


•• ■€>—■: 




















\ 




9 




•••-.©• 3 













0.4 0.6 

r eff 



(a) 



0.0 
-1.0 
^ -2.0 
^ -3.0 
1 -4.0 
uf -5.0 
-6.0 
-7.0 









w/o $ 
with $ 








' "C 












>. 













0.4 0.6 

r eff 



(b) 







w/o 
with 










\ 

























0.2 0.4 0.6 0.E 
r eff 



(c) 



Figure 3: Particle flux of C 6+ (a), Ne 8+ (b) and Fe 20+ (c) as a function of r e ff in LHD 
for the set of profiles considered in case I (see text). 
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Figure 4: Particle flux of C 6+ (a), Ne 8+ (b) and Fe 20+ (c) as a function of r e ff in LHD 
for the set of profiles considered in case II (see text). 



fig. [2](b) the $ map dependent on r e ff and 9 for case I is given, while [2](c) represents the 
one corresponding to the case II. 

The simulations run by EUTERPE were carried out for C 6+ , Ne 8+ and Fe 20+ . The radial 
profiles of the particle flux are shown in figs. [3(a)- (c) for case I. In this set of plots 
the impact of $ on the three species is detrimental up to r e ff ^0.7 from the impurity 
accumulation perspective. From that radial location outwards the trend is reversed and 
the inward flux in the presence of $ becomes weaker. We conclude that $ can act either 
to amplify or mitigate for the inward flux driven by E r . 

The same behaviour is shown for case II in figs. ffl(a)-(c). Again l> breaks the monotonic 
growth of the inward impurity flux with increasing r e ff. The change in the trend is 
observed at a different radial location than in case I, r e ff ~ 0.4, which corresponds 
to almost the position where the electrostatic potential starts to be appreciably large, 
see figure [2(c). A remarkable feature in this case is that the particle flux is zero at 
r e ff ~ 0.7 for C 6+ and Ne 8+ and considerably less negative for Fe 20+ . That radial 
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Figure 5: (Left) Fourier coefficients $^ and $^ corresponding to the $ map shown in 
plot [2](b). (Right) Fourier coefficients of $ corresponding to the plot[2](c). In Fourier 
space $ is in GSRAKE expressed as <£> = ($^ sin m9 + $^ cos mO) with m the poloidal 
mode number. 



position corresponds to the maximum amplitude of $ along 9. 

At this point it is convenient to recall the discussion about the role of $ for driving 
radial transport and trapping particles. First, it is important to notice that in case 
I the variation of $ leads to an electrostatic energy variation of approximately 1% of 
the thermal energy for a unit charge. In case II the same variation would result in 
approximately a 2% change. This indicates that the strong change in the behaviour of 
the particle flux is not due to the increase of the ratio q<&/T only. Looking at eq. ([3j) 
and eq. (Ill) we can conclude, however, that for typical values of AB/B and R/a in the 
configuration under study v$ and a$ can be comparable to Vd and a m respectively. This 
means that $ can drive as much impurity transport as the inhomogeneity and curvature 
of the magnetic field, and also can trap as many impurities as the magnetic helical wells. 
Thus, the relative position of the magnetic field and electrostatic wells may counteract 
the inward flux driven by E r . 

Indeed, if we look at the Fourier coefficients of $ in fig. [5] it can be observed that the 
mitigation of the inward flux coincides in both cases I and II with the radial location 
where the sin(0) component is relatively low, at the edge in case I and across the plasma 
in case II. 

4. The self-consistent calculation of $ (0, 0, t) in EUTERPE 

In light of these results it is of immediate interest to extend the calculations to account 
for the toroidal dependence of $. Apart from the role that this can play, it is important 
to remember that the ripple average also reduces the number of possible configurations 
than can be studied. Although this is beyond the scope of the present work, we show 
in this section a preliminary result where 0, t) is obtained self-consistently with 
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EUTERPE. The approach rests on the fulfilment of the neutrality condition on the density 
perturbations for each species, J2 S ^ = 0- The starting point is the quasi-neutrality 
gyro- kinetic equation [12] that the code solves at each time step. Considering for 
simplicity kinetic ions and adiabatic electrons, the equation reads as follows: 

e (fk) + m, V • (|| VjS) - ^ ($-$)= 0, (12) 

where the lower indices i and e denote bulk ions and electrons respectively, n is the 
equilibrium density, (n^) is the gyro- averaged ion density, $ is the flux-surface averaged 
electrostatic potential, and rrii the ion mass. The equation is simplified by invoking the 
limit k±p <C 1 with k± the perpendicular characteristic wave-length of the fluctuations. 
Since the second term in eq. ( IT2l ) is a factor p 2 k\ smaller than the others it can be 
neglected. In addition, the gyro-average operation is unnecessary in drift kinetics. And 
finally, by flux surface averaging of the remaining expression, it can be shown that 
$ = $o,o = 0. The sub-indices in $ denote that the poloidal and toroidal mode numbers, 
m and n respectively, are both equal to zero. The final expression is then: 

6 

hi = — n 0e $. (13) 

J- e 

In order to prevent the code from developing short wave-length unstable modes, only 
low mode numbers are retained by employing a spectral filter set to < m < 4 and 
-4 < n < 4. 

Since the trajectories are radially local, it is difficult to calculate E r . An iterative 
adjustment of E r until ambipolarity is fulfilled is not feasible due to the dynamic nature 
of the simulation, which makes it impossible to determine the flux before it has saturated. 
Thus, E r is imposed externally. 

Figures [6](a)-(b) show the result of a simulation for !>(#,(/>, i) at r e ff = 0.5, considering 
the profiles of the labelled as case II in section El The calculation of both the hydrogen 
flux Fi and $ is dynamical, as has already been mentioned. In figure [6](a) Fi is shown 
as a function of time. Once it has reached the stationary value, a time average of the 
potential is performed. The resulting time averaged potential is represented 

in fig. [6](b). The time interval considered for the average is coloured yellow in fig. [6](a). 
It can be noticed that a weak but appreciable variation of $ with (j) is present. 

5. Summary and discussion 

In the present work we have tackled the problem of impurity particle transport in 
stellarator geometry. The motivation is that impurity accumulation can occur under 
certain plasma conditions and is predicted by the standard local and mono-energetic 
neoclassical theory. 

Nevertheless, the mono-energetic approximation does not hold for impurities, whose 
high charge makes them sensitive to electrostatic potential variations within the flux 
surface. This makes the problem at least - assuming radial locality of the trajectories in 
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Figure 6: (Left) Time dependent flux of hydrogen in LHD for the case II profiles (see 
text in section [3]) at r e ff = 0.5. (Right) Time averaged electrostatic potential 0). 

the lowest order - a 4D problem. The task of solving it for the impurity particle flux F z 
has been considered by the Monte-Carlo code EUTERPE. The potential $(#) arising by 
enforcing quasi-neutrality in the code GSRAKE was used as input for the EUTERPE runs. 
The results have shown that $ can both increase and decrease the impurity 
accumulation, depending on the interplay between the radial drives v$ and Vd and 
the electrostatic and magnetic trapping. These, in turn, have shown to be determined 
ultimately by the Fourier spectrum of $ rather than by its absolute value. The 
calculations show that the inward flux of impurities can be suppressed completely. 
Future work aims at a self-consistent calculation of <E> by EUTERPE along the line that 
has been presented. The preparation of a global neoclassical version of the code, which 
would allow E r to be computed, is also forthcoming. 
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